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ABSTRACT 

When one considers the effect in the physical space, Daubechies-based wavelet 
methods are equivalent to finite difference methods with grid refinement in regions of 
the domain where small scale structure exists. Adding a wavelet basis function at a 
given scale and location where one has a correspondingly large wavelet coefficient is, 
essentially, equivalent to adding a grid point, or two, at the same location and at a 
grid density which corresponds to the wavelet scale. This paper introduces a wavelet- 
optimized finite difference method which is equivalent to a wavelet method in its 
multiresolution approach but which does not suffer from difficulties with nonlinear 
terms and boundary conditions, since all calculations are done in the physical space. 
With this method one can obtain an arbitrarily good approximation to a conservative 
difference method for solving nonlinear conservation laws. 
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1 Introduction 


In the numerical simulation of equations which model physics it is common that small 
scale structure will exist in only a small percentage of the domain. If one chooses 
a uniform numerical grid fine enough to resolve the small scale features then in the 
majority of the domain the solution to the equations will be over resolved. One would 
like, ideally, to have a dense grid where small scale structure exists and a sparse grid 
where the solution is composed only of large scale features. 

Consider now a Daubechies-based wavelet system. Wavelets provide a natural 
mechanism for decomposing a solution into a set of coefficients which depend on scale 
and location. One can then work with the solution in a compressed form where one 
works only with the wavelet coefficients which are larger in magnitude than a given 
threshold. Wavelets, therefore, sound ideal for solving the type of problem mentioned 
the previous paragraph. There are, however, serious problems matching the order 
of differentiation accuracy at the boundary for nonperiodic boundary conditions, see 
[9], with the superconvergence encountered with periodic boundary conditions, see 
[7]. Furthermore, wavelet methods generally require a tranformation between the 
physical space and the coefficient space for either evaluation of nonlinear terms or for 
differentiation. 

In this paper a wavelet method which satisfies the goals of the first paragraph 
while using the wavelet machinery outlined in the second paragraph without the 
accompanying difficulties encountered at the boundaries and the expense of constantly 
tranforming between the physical space and the coefficient space will be introduced. 
That is, the new method utilizes the strength of wavelets, scale detection and data 
compression, while avoiding the difficulties by using wavelets in their finite difference 
form. 

The following is a list of the sections of this paper with the noteworthy points. 
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1. Introduction 


2. Wavelet Definitions and Notation 

3. Finite Difference Grid Refinement and Wavelets: 

This section will establish that Daubechies-based wavelet methods are equal to 
finite-difference methods with grid refinement. 

4. The Wavelet-Optimized, Adaptive Grid, Finite Difference Method: 

In this section a new numerical method which utilizes the strength of wavelets 
and avoids the difficulties will be proposed. That is, wavelets will be utilized 
to define the grid for finite difference methods. The new method is named, 
‘The Wavelet-Optimized, Adaptive Grid, Finite Difference Method’, or simply 
‘WOFD’. 

5. WOFD applied to Burgers’ equation: 

Numerical results of the wavelet-optimized, adaptive grid, finite difference method 
applied to Burgers’ equation with periodic and nonperiodic boundary conditions 
will be given. 

6. Accuracy of WOFD: 

The error in finite-difference derivative approximation on a 5-point stencil is of 
the form, 

Err = AiA 2 A 3 A 4 ^/ M (a). 

Think of f(x) as a pure mode, f(x) = e lx ^, where £ is frequency or wave number. 
When the data is locally smooth, i.e., composed of low frequencies, the wavelet 
coefficients are small and consequently the A’s are allowed to be large. When 
the data is locally oscillatory, i.e., composed of high frequencies, the wavelet 
coefficients are large and WOFD reduces the size of the A’s. The effect is 
that the derivative approximation error will not grow faster than linearly with 
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respect to frequency. Recall that without grid adjustment this error would grow 
as a fifth power of frequency for a fourth order scheme. 

7. Stability of WOFD: 

Analytical stability methods are beyond reach due the arbitrary nature of the 
grid. But, in practice the method displays no instability when applied to Burg- 
ers’ equation. 

8 . Efficiency of WOFD: 

The WOFD method finds an approximation to the solution found on the finest 
scale across the whole domain. The efficiency depends on the rate of data 
compression. That is, if the finest scale has N grid points and the WOFD 
averages, say, N 0 grid points, then the WOFD method will fold the solution 
using, roughly, ^ times the number of operations used to fold the finest grid 
solution. 

9. WOFD in Higher Dimensions 

The discussion here will be limited to grid selection, ft will be seen from a few 
examples that WOFD is an effective method for grid selection in higher dimen- 
sions. The examples given are for two dimensions. The local ‘spectral analysis’ 
of a wavelet method provides exactly the information needed to thoroughly 
understand the data and, hence, define a grid properly. 

10. Conclusion: 

The WOFD method is an efficient and stable alternative to a Daubechies wavelet 
method. The WOFD method and a wavelet method are essentially the same. 
The only significant difference is the manner in which the grid is refined. The 
WOFD method, by contrast, avoids difficulties with nonlinear terms and bound- 
aries by performing all calculations in the physical space. 


3 



2 Wavelet Definitions and Relations 


The term wavelet is used to describe a spatially localized function. ‘Localized’ means 
that the wavelet has compact support or that the wavelet almost has compact sup- 
port in the sense that outside of some interval the amplitude of the wavelet decays 
exponentially. We will consider only wavelets that have compact support and that 
are of the type defined by Daubechies [4] which are supported on [0,2 M — 1], where 
M is the number of vanishing moments defined later in this section. 

To define Daubechies wavelets, consider the two functions <f>(x) and ip(x) which 
are solutions to the following equations: 


where <f>(x) is normalized, 


Let, 


and 


L—l 


(j>{x) = \[2 ^2 h k (f)(2x - k), 

k= 0 
L—l 

(1) 

= \/2 gk<P{2x - k), 

k= 0 

/*oo 

(2) 

/ ( t>(x)dx = 1. 

J — CO 

(3) 

<f>{(x) = 2 - 2 (f){2~ 3 x — k), 

(4) 

■ ipl(x ) = 2 ~ 2 il>( 2 ~ 3 x — k), 

(5) 


where j,k £ Z , denote the dilations and translations of the scaling function and the 
wavelet . 

The coefficients H = { ^ ^ = o an d bf = {dk}k = o are related by g k = ( — 1 ) k liL_k for 
k = 0, L — 1. Furthermore, H and G are chosen so that dilations and translations 
of the wavelet, if) 3 k (x) } form an orthonormal basis of L 2 (R ) and so that ip(x) has M 
vanishing moments. In other words, ip 3 k (x ) will satisfy 

/ OO 

( 6 ) 

-OO 
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where Ski is the Kronecker delta function. Also, ip(x) = ifio( x ) satisfies 



for m = — 1. Under the conditions of the previous two equations, for any 

function f(x) £ L 2 (R ) there exists a set {djk} such that 


f( x ) = djki’Kx), (s) 

jez Kez 

where 

/ OO 

f( x Wk( x ) dx - ( 9 ) 

-oo 

The number of vanishing moments of the wavelet ip(x) dehnes the accuracy of 
approximation. The two sets of coefficients H and G are known in signal processing 
literature as quadrature mirror filters [5]. For Daubechies wavelets the number of 
coefficients in H and G, or the length of the filters H and G, denoted by L, is related 
to the number of vanishing moments M by 2 M = L. For example, the famous Haar 
wavelet is found by defining H as h 0 = hi = 1. For this filter, H, the solution to 
the dilation equation (1), <f>(x), is the box function: <f>(x) = 1 for x £ [0,1] and 
<f>(x) = 0 otherwise. The Haar function is very useful as a learning tool, but it 
is not very useful as a basis function on which to expand another function for the 
important reason that it is not differentiable. The coefficients, H , needed to define 
compactly supported wavelets with a higher degree of regularity can be found in [4], 
As is expected, the regularity increases with the support of the wavelet. The usual 
notation to denote a Daubechies wavelet defined by coefficients H of length L is Dl. 

It is usual to let the spaces spanned by <f>k( x ) and V’fclu) over the parameter k, 
with j fixed, to be denoted by Vj and Wj respectively: 


t t span ,j , . 

Vj = 9k\ x ) 

J kez y 
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The spaces Vj and Wj are related by [4] 

... C b C ho C V-i C ..., (12) 

and that 

r, = ii +1 ®in +1 . (13) 

The previously stated condition that the wavelets form an orthonormal basis of L 2 (R ) 

can now be written as, 

L 2 (i?) = ©»'. (14) 

iez 

Two final properties of the spaces Vj are that 

n ^ = {»}■ (1 5 ) 

j ez 

and 

\JV S =L\R). (16) 

Of course, inhnite sums and unions are meaningless when one begins to implement 
a wavelet expansion on a computer. In some way one must limit the range of the 
scale parameter j and the location parameter k. Consider first the scale parameter j . 
As stated above, the wavelet expansion is complete: L 2 (R ) = ©jgz^j- Therefore, 
any f(x) £ L 2 (R ) can be written as, 

f( x ) = d l^l( x ), 

jez kez 

where due to orthonormality of the wavelets d J k = f(x)if) 3 k (x). In this expan- 

sion, functions with arbitrarily small-scale structures can be represented. In practice, 
however, there is a limit to how small the smallest structure can be. This would 
depend, for example, on how fine the grid is in a numerical computation scenario or 
perhaps what the sampling frequency is in a signal processing scenario. Therefore, 
on a computer an expansion would take place in a space such as 

y 0 = w k © W 2 © . . . © Wj © Vj, (17) 
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and would appear as, 


p vj{x) = Y, s k^k( x ) + YJ2 d U’k( x )> ( 18 ) 

kez j= 1 kez 

where again due to orthonormality of the basis functions d k = and 

s k = f(x)cj) J k (x). In this expansion, scale j = 0 is arbitrarily chosen as the finest 
scale that is needed, and scale J would be the scale at which a kind of local average, 
cf) k (x), provides sufficient large scale information. In language that is likely to appeal 
to the electrical engineer it can be said that 4> J k {x) represents the direct current portion 
of a signal at location k and that ip k ( x ) represents the alternating current portion of 
a signal at, very roughly, frequency j and location k. As stated above, one must 
also limit the range of the location parameter k If one assumes periodicity, then the 
periodicity of f(x) induces periodicity on all wavelet coefficients, s k and d k , with 
respect to k. Without periodicity, the location parameter k will begin at 1 with the 
left-hand side boundary functions and end with some maximum number N at the 
right-hand side boundary functions. 

This completes the definition of wavelets. 
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3 Finite Difference Grid Refinement and Wavelets 


In this section it will be shown that Daubechies-based wavelet methods when con- 
sidered in the physical space are equivalent to explicit finite difference methods with 
grid refinement. In a Daubechies wavelet method the ‘refinement’ is accomplished by 
adding wavelet bases functions in regions where structure exists corresponding to the 
scale of the wavelet used for analysis. In a finite difference method the ‘refinement’ is 
accomplished by adding grid points in regions chosen by some grid refinement mech- 
anism. In this section it is argued that since wavelet methods correspond to central 
finite difference operators when the grid is uniform and since wavelet methods contain 
a natural and effortless mechanism for ‘grid refinement’, then one can simply use the 
wavelets to refine a grid for finite difference operators. In this way one can maintain 
the superconvergence encountered with periodic boundary conditions, see [7], which 
is lost when one constructs wavelets on an interval, see [9]. That is, boundary condi- 
tions are imposed in the same manner as for finite difference operators. Furthermore, 
there is no longer a difficulty with nonlinear terms requiring constant transformation 
between the physical space and the coefficient space since all calculations are done in 
the physical space. 

This section contains four subsections: 

1. The wavelet decomposition matrix will be constructed. 

2. It will be seen that under the assumption of periodicity and without data com- 
pression that the effect in the physical space of differentiation in the D 4 subspace 
Vo is exactly the same as differentiation with the optimal central 4th-order finite 
difference operator. 

3. Now we compare wavelets and finite difference in the subspace Vo = W\ ® V\. If 
Ax is the grid spacing in Vo then 2Ax is the grid spacing in V\ and the wavelet 
coefficients in W\ indicate if refinement is needed for a local grid spacing of Ax. 
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4. Finally, a division of the subspace Vo which might be used in practice is studied: 
Vo = Wi ® W 2 ® W 3 ® V3. Similar to above, the grid spacing in the subspace V3 
would be 8 Ax. The first refinement is controlled by the subspace W 3 which can 
refine locally to a grid spacing 4 Ax. Likewise, the subspace W 2 refines locally 
to a grid spacing of 2 Ax and W\ to a local grid spacing of Ax. 

3.1 Wavelet Decomposition Matrix 


The wavelet decomposition matrix is the matrix embodiment of the dilation equation 
defining the scaling function and the accompanying equation defining the wavelet. 
The following two recursion relations for the coefficients s k and d J k can be found from 
equations (1) and (2), respectively: 

n=2M 

4 = I] 44+4-2, ( 19 ) 

n=l 


and 


n=2M 

4 = I] ^4+4-2' (2°) 

n=l 

Denote the decomposition matrix embodied by these two equations, assuming peri- 
odicity, by 44xiv w here the matrix subscripts denote the size of the matrix, and the 
superscripts indicate that P is decomposing from scaling function coefficients at scale 
j to scaling function and wavelet function coefficients at scale j + 1. Let Sj contain 
the scaling function coefficients at scale j. (Note that when vector notation is used 
the scale is given as a subscript.) P therefore maps Sj onto and dj + \\ 


pjj+i 

r NxN 


®+l 

4+i 


( 21 ) 


Note that the vectors at scale j + 1 are half as long as the vectors as scale j. To 
illustrate further, suppose the wavelet being used is the four coefficient D 4 wavelet, 
and suppose one wants to project from 8 scaling function coefficients at scale j to 4 
scaling function coefficients at scale j + 1 and 4 wavelet coefficients at scale j + 1. 
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The decomposition matrix in this case is, 


Pi 


j'd+i — 


8x8 — 


hi h2 h 3 hi 
0 0 hi h 2 

0 0 0 0 
h 3 hi 0 0 

9i 92 93 9i 
0 0 gi g 2 g 3 94 

0 0 0 0 
J3 «4 0 0 


0 

h 3 

hi 

0 

0 


0 

hi 

h 2 

0 

0 


0 

0 

h 3 

hi 

0 

0 


0 

0 

hi 

h 2 

0 

0 


9 1 9 2 93 94 

0 0 ft j 2 


( 22 ) 


where the periodicity is seen from the coefficients ‘wrapping around’. 

Now let us consider differentiation. Let the four matrices A 3 NxN} B 3 NxN} Cj [ 


NxNi 


and R 3 NxN} see [7] and [1], contain the derivative projection coefficients, 


A J : dj 


d 


J 5 


V • 


dj i 


C 3 : di — > s 


‘i 


R 3 : s 


J 


where Sj and dj denote the coefficients of the expansion of the derivative of a function 
which is initially defined by the expansion coefficients Sj and dj. The exact form of the 
matrices A, B } and C is not important for the discussion here. The important point is 
the form of the matrix R. It is always a finite difference operator. For the Di wavelet 
R corresponds to the optimal central 4th-order finite difference operator. For higher 
order wavelets, Z) 6 , D 8} etc., R is a finite difference operator, but it is not optimal 
in the sense of using the minimum number of coefficients to obtain a given accuracy. 
The numerical values of the coefficients were found in [1] and the superconvergence 
accuracy was proven in general in [7]. For the Di wavelet an explicit 8x8 example 
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of matrix R is, 


-RsyS — 


_ 2 
3 

x 

12 

0 


12 

2 

3 


2 

3 

0 

_ 2 
3 

x 

12 

0 

0 

0 

_x 

12 


12 

2 

3 

0 

_ 2 
3 

x 

12 

0 

0 

0 


12 

2 

3 

0 

_ 2 
3 

x 

12 

0 

0 


12 

2 

3 

0 

_ 2 
3 

x 

12 

0 


12 

2 

3 

0 

_ 2 
3 

x 

12 


x 

12 

0 


0 
0 

0 0 

--T 0 


12 

2 

3 

0 

_ 2 
3 


_ 2 
3 

x 

12 

0 

0 

0 


12 

2 

3 

0 


( 23 ) 


We will now see how ‘grid refinement’ is accomplished in a wavelet scenario by 
examining three divisions of the subspace Vo in three following three subsections. 


3.2 Wavelet Expansion and Derivative in Vq 


One can calculate the derivative of a wavelet expansion at any level in the wavelet 
decomposition. This subsection will explore the first of three of the options. To be 
explicit, suppose that a periodic function f(x) has been approximated on a grid with 
16 scaling function coefficients to get s 0 , and for the current argument assume that 
the coefficients have been calculated exactly. Note that periodicity of f(x) induces 
periodicity on the coefficients s 0 . The coefficients of the expansion of yy/(®) in Vo 
are found from s 0 by an application of the matrix -Ri 6x16 : 


1 

0 i—i 
<0 

1 


1 

0 i—i 
v t0 

1 

s ° 

*2 


s ° 

*2 

S° 

*3 


S° 

*3 

S° 


i° 

*4 

S° 


i° 

*5 

S° 


S° 

*6 

S° 

*7 


4° 

*7 

s ° 

1 pO 

Ax ^16 X16 

4° 

*8 

S° 


i° 

*9 

S° 

*10 


i° 

*10 

s ° 

*11 


4° 

*n 

s ° 

*12 


4° 

*12 

S° 

*13 


4° 

*13 

S° 

*14 


4° 

*14 

S° 

*15 


i° 

*15 

1 

to 

O i—i 
<0 

1 


1 

to 

O i—i 
v t0 

1 


( 24 ) 
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Let us now examine the entire process of going from point values in the physical 
space to scaling function coefficients in Vo, differentiating, and finally returning to 
point values of the differentiated function in the physical space. Begin with f(x) £ 
L 2 (R) defined at 16 evenly-spaced points over [— 7r,7r) and let f(x) is 2tt periodic. 
To differentiate the samples of f(x), /, with the 4-th order optimal central finite 
difference operator, say DfdA, we get, 

/ = D m f. (25) 

Now, suppose that we have mapped these 16 samples into the scaling function coeffi- 
cients in Vo by applying the circular, periodicity implies circularity, see [7], quadrature 
matrix Q, 

■So = Qwxwf- (26) 

We now fold i 0 by applying ^yi?i6xi6- The two matrices R and Q are, however, both 
circular and, hence, commute: 

■So = Qi6xi6~: Rwxwf- (27) 

I\x 

Now, returning to the physical space we get, 

/ = Qi6xi6Qiexi6~r—Ri6xi6f } (28) 

I\x 

and we are back to equation (25) again since, 

Dffa = — iWe- (29) 

Hence, we have shown that under the assumption of periodicity and without 
data compression that the _D 4 wavelet differentiation corresponds exactly to optimal 
central 4th-order finite differencing. Note that data compression is the goal of any 
wavelet method. The embodiment of data compression in the physical space is a 
nonuniform grid. That is, the grid must be dense in regions where small structure 
requires foie resolution and the grid can be sparse when the data is composed of large 
scale components. 


12 



Now we move to the first decomposition of Vo = W\ ® Vi in which data compression 


can be achieved. 

3.3 Wavelet Expansion and Derivative in W\ © V\ 


Consider now a decomposition of the vector of scaling function coefficients s 0 onto 
the scaling function and wavelet coefficients at scale j = 1 by an application of the 
matrix P^xie- 


O i—i 
<0 


’ 4 ’ 

s° 


4 

S° 

*3 


s 1 

*3 

s° 


s 1 

*4 

S° 


A 

*5 

S° 


4 

S° 

*7 


4 

s° 

p0,l 
^16 X 16 

s 1 

L *8 J 

S° 


d\ 

S° 

*10 


4 

s° 

*11 


Pi- 
co h-i 

s° 

*12 


d\ 

S° 

*13 


d\ 

S° 

*14 


4 

S° 

*15 


d\ 

to 

O i—i 
<0 


CO h-i 


As in Vo, we have 16 basis functions, but now the subspace Vo is decomposed into 


‘low frequency’, Vi, and ‘high frequency’, IV© components: Vo = Vi ® IV© In order 
to calculate the coefficients of the derivative expansion in Vi ® W\ the following 
projections are calculated: 


■Sl 

d\ 


1 

Dl 

-^8x8 

fi 

^8x8 


© 

2Ax 

pi 

-°8x8 

A 1 

^*-8x8 


di 


(31) 


If one now applies the matrix (C’ 1 °g! <16 ) T (T denotes transpose and hence inverse for 
this unitary matrix) to the derivative coefficients at scale j = 1 one gets, 


[ 4 ] = (ifixief ' 

and one gets exactly the same coefficients as before when the matrix ^r-Ri6 Xl 6 
applied to s 0 . 



(32) 


was 
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Now suppose that /(x) is smooth enough such that a grid of eight points provides 
sufficient resolution. Define / 2 to be the 8 element vector containing every other entry 
of the 16 element vector /. 4-th order differentiation of / 2 is performed by applying 


2Ax^ 8x8 ’ 


/ 2 = 


1 

2A: 


'©8x8/2- 


(33) 


Similar to above, we project the eight dimensional / 2 into the eight dimensional 
wavelet subspace Vi using Q 8x s and differentiate to get, 


© 


2A: 


-R%xsQsxsf2i 


(34) 


followed by projection back into the physical space with the matrix Qg^s we get 
equation (33) again. 

That is, we have seen that if we work only in Vo that we have 4th-order finite 
differencing with a grid spacing of Ax, whereas if we work only in Vi we have 4th- 
order finite differencing with a grid spacing of 2 Ax. But, the two subspaces Vo and Vi 
are related by Vo = Vi ® Vfi. Recall, that the subspace Vfi contains bases functions 
which are locally oscillatory and are compactly supported. An inner product of these 
bases with the data f(x) will detect local oscillations in f(x) and provide exactly the 
information necessary to refine the grid locally from 2 Ax to Ax. This wavelet grid 
refinement mechanism can be used not only to add wavelet bases functions where one 
has a large inner product but also to refine the grid in the same region and at a scale 
corresponding to the wavelet scale. 


3.4 Wavelet Expansion and Derivative in W\ © W 2 © W 3 © V 3 

Let us close this section with a wavelet decomposition that one might use in practice. 
That is, again Vo denotes the finest scale subspace and we decompose Vo as, 


Vq — W\ ® IV2 ® W3 ® V3. 


(35) 
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The vector of coefficients in this subspace is obtained by the following decompositions: 



Let us suppose that we have performed the above wavelet decomposition on a 
vector of data points, /, at some point during a simulation which contains data 
at many different scales. Furthermore, let there be a shock, or a near shock, near 
the right-hand boundary. The coefficients and s \ represent local averages in the 
subspace V 3 corresponding to the ‘base grid’ of size 8 Ax and will not yield much useful 
information with respect to the shock. The coefficients d\ and d \ of the subspace W3 
will yield the presence of oscillations of relatively large scale. A true shock contains 
all frequencies and one would expect to have some coefficient perturbation even in 
W 3 yielding grid refinement to a grid spacing of 4Ax in a neighborhood of the shock. 
The coefficients d 2 , d%, d 2 } and d 2 A in the subspace LF 2 will detect oscillations at the 
corresponding scale only near the shock. That is, the first two coefficients d 2 and d\ 
are responsible for detecting small scale structure at the left-hand side of the domain 
which is away from the shock, and we, therefore, expect that they will be near zero 
in magnitude. The coefficients d 2 and on the other hand, are positioned near 
the shock and will have a relatively large amplitude indicating the presence of small 
scales. The grid will, therefore, be refined to a spacing of 2Ax at the right-hand side 
of the domain. Likewise, the remaining coefficients in the subspace W\ will refine the 
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grid to a spacing of Ax at the right-hand side of the domain. 

In conclusion, this section has been devoted to first illustrating how Daubechies- 
based wavelet methods are in essence finite difference methods with grid refinement, 
and second to illustrating how the Daubechies-based wavelets can be used to de- 
fine a grid for finite difference methods. The next section will make this symbiotic 
relationship between Daubechies wavelets and finite difference methods formal. 
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4 The Wavelet-Optimized, Adaptive Grid Finite 
Difference Method, (WOFD) 

The new method is to apply finite difference on a grid which is defined by the magni- 
tude of wavelet coefficients at various scales. That is, wavelets can detect oscillations 
in a function at any location and scale. Given a function f(x) for x £ /, where / is 
some interval, one decomposes f(x) into a set of wavelet coefficients which depend 
on two parameters, one for location and one for scale, say d J k , where k is the location 
parameter and j is the scale parameter. If a wavelet coefficient is large in magnitude, 

I4l > T. (37) 

or large in energy (In practice the two criteria yield roughly the same grid.), 

«) 2 > T, (38) 

where T is a coefficient threshold chosen by the user, then WOFD adds a grid point, 
or two, at location k and at a grid density corresponding to the scale j. That is, 
WOFD defines a grid which will completely resolve a function across the entire domain 
without over resolving it where it is relatively smooth, or composed only of large scale 
structure. For the specific case of the D 4 wavelet outlined in the previous section, the 
D 4 wavelet decomposition provides the optimal grid for 4th-order finite differencing. 

The grid definition should be made by a Daubechies wavelet which corresponds in 
terms of superconvergence accuracy to the accuracy of the finite difference operator. 
That is, it was proven in [7] that the differentiation matrix for the Daubechies wavelet 
D 2 M 1 where M is the number of vanishing moments, displays differentiation accuracy 
of order 2 M under the assumptions of periodicity and a uniform grid. Recall, that this 
wavelet subspace can only represent exactly the first M polynomials as determined 
by the number of vanishing moments. This order of accuracy 2 M should equal the 
order of accuracy of the finite difference operator for optimal grid selection. 

In the next section WOFD will be applied to Burgers’ equation. 
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5 WOFD Applied to Burgers’ Equation 


In this section WOFD will be applied to Burgers’ equation, 

dU_ _ dU_ d 2 U 

dt dx dx 2 ’ 


with the initial condition, 


u(x, 0) 


x Lu , , 

— | — simz-KX 

3 3 V 


( 39 ) 


(40) 


The goal of this section is to illustrate that WOFD using the D 4 wavelet produces 
a solution on a nonuniform reduced grid which is ‘equivalent in character’ to the 
solution provided by 4th order finite differencing on the finest uniform grid. That 
is, for a given viscosity, e, there exists a grid size fine enough such that oscillations 
do not develop at the ‘shock’. This can be made more precise by saying that one 
has a grid fine enough such that the total variation of the solution does not increase. 
‘Equivalent in character’ means that the total variation of the solution provided by 
WOFD increases if and only if the total variation of the solution produced by finite 
differencing on the finest uniform grid increases. 

In all the following plots the uniform finite differencing is provided by the optimal 
central 4th-order finite difference operator. The temporal discretization is achieved 
by 4th-order Runge-Kutta. The WOFD coefficient threshold which determines which 
grid points to use based on the wavelet coefficient magnitude is set to T = .001. 
Note that when the WOFD coefficient threshold is set to T = 0 that one gets finite 
differencing on the uniform finest grid. In addition, if the coefficient threshold is set 
to a very large number, say T = 100, then one gets finite differencing on a uniform 
sparse grid. The size of this sparse grid is determined by the number of wavelet 
decompositions one specifies. 


5.1 Periodic Boundary Conditions 

In figure (1) WOFD is compared to finite differencing on uniform grid sizes 32, 64, 
and 128. The upper left-hand plot has the WOFD solution superimposed on the finite 
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difference solution. The two solutions are visually indistinguishable. Along the x-axis 
of the plot an ‘x’ is placed at every position where a WOFD grid point is used. One 
can see that the grid points are dense at the shock and sparse where the solution is 
smooth. 

The remaining three plots show finite difference solutions on various uniform grid 
sizes. One sees oscillations for grid sizes 32 and 64 but not for grid size 128. 

Figure (4) provides an additional plot for periodic boundary conditions with a 
slightly larger viscosity. 

5.2 Nonperiodic Boundary Conditions 

The boundary conditions considered here are such that the boundary values of the 
solution are required to be fixed at the initial condition values, 

u(0 } t) = u(l,t) = i. (41) 

In all the plots, differentiation at the boundary for the uniform finite difference method 
is achieved by the optimal 4th order one sided finite difference coefficients, see [2], 
For WOFD both 4th order and 3rd order boundary differentiation will be considered. 

In Figure (2) the differentiation at the boundary is 4th order, and, as above, 
WOFD provides a solution which is ‘equivalent in character’ to the finite difference 
solution on the finest grid, 128 grid points, while reducing the number of degrees-of- 
freedom necessary to achieve this solution. 

In Figure (3) the differentiation at the boundary is 3rd order. The solution for the 
3rd-order boundary differentiation is good, but a slight difference can be seen with 
the finite difference method on the finest grid. Again, finite difference on more coarse 
grids oscillates more at the shock than the finest grid solution. 

Figure (5) provides an additional plot illustrating the solution provided by WOFD 
for the nonperiodic case. 
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FDon Grid 128 and WOFD 


FD on Grid 128 



Figure 1: Illustration of equivalence of WOFD to an equivalent order finite difference 
method applied across the entire domain at the finest scale. The boundary conditions 
are periodic. Final time = 2, Viscosity = .02. 
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FDon Grid 128 and WOFD 


FD on Grid 128 




FD on Grid 64 FD on Grid 32 




Figure 2: Illustration of .equivalence of WOFD to an equivalent order finite difference 
method applied across the entire domain at the finest scale. The boundary values 
are fixed at the initial condition values. Differentiation at the boundary is 4th order. 
Final time = .3, Viscosity = .005. 


21 







FDon Grid 128 and WOFD 


FD on Grid 128 




FD on Grid 64 FD on Grid 32 




Figure 3: Illustration of .equivalence of WOFD to an equivalent order finite difference 
method applied across the entire domain at the finest scale. The boundary values 
are fixed at the initial condition values. Differentiation at the boundary is 3rd order. 
Final time = .3, Viscosity = .005. 
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6 Accuracy of WOFD 


In this section the order of accuracy will be examined. For numerical methods where 
the grid is uniform the order of accuracy is clearly defined. For WOFD, on the other 
hand, the discussion of accuracy is slightly more complicated. That is, it can be said 
that WOFD approximates a 4th-order finite difference method as well as one desires, 
and when the coefficient threshold is set to zero then WOFD is truly 4th-order. So, 
WOFD approximates methods of a given order as well as is desired. In addition, it 
will be seen that the WOFD method has a very nice feature that the rate of growth 
of the error in approximating the derivative is at most a linear function of frequency. 

6.1 Error in Derivative Approximation 

The finite difference equations used in this paper use either a 4-point stencil or a 
5-point stencil. The derivative approximation error for the equations on a 3-point 
stencil will be given. The derivative approximation error for a larger stencil is an 
obvious extension of the error given here. 

Consider the following Lagrangian interpolation of a quadratic polynomial through 
the three points: (aq, /(aq)), (x 2 , f(x 2 )), and (x 3 , f(x 3 )), for x 1 < x 2 < x 3 : 

g(x) = (42) 

(x - x 2 )(x - x 3 ) (x - x-l)(x - x 3 ) (x - x-l)(x - x 2 ) 

1 (aq - X 2 ){x x - x 3 ) 2 (x 2 - x 1 )(x 2 - x 3 ) 3 (x 3 - aq)(aq - x 2 ) 

If we differentiate g(x) and evaluate at x 2 we get, 

^ g ( x ) U = + A 1 A 2 i/("')(a), (43) 

for some a £ [aq, aq], where Ai = x 2 — x 1 and A 2 = x 3 — x 2 . 

6.2 Control of Error Growth 

As given above we will examine the special case of a 3-point stencil where the grid is 
refined by the Haar wavelet. In practice I never use the Haar wavelet, but it is very 
useful as an illustration tool. 
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As given above the error is, 


Err = AiA 2 — H a )- (44) 

for some a £ [aq, x 3 \. Now, let / be a pure sinusoid of frequency £: f(x) = e*W Then 
the magnitude of the error becomes, 

6 1 -Err | = AiA 2 ^ 3 . (45) 

The grid refinement mechanism used by WOFD is such that, roughly, 

A = 1 (46) 

The magnitude of the error becomes, 

6|Err| = £. (47) 

That is, the refinement mechanism keeps the rate of growth of the error linear with 
respect to frequency. Whereas, without the grid refinement the error grows as a cubic 
in this case. 

This is one particular refinement mechanism, but is representative of a typical 
refinement method. 

6.3 Relationship of Threshold Size to Solution 

The grid for WOFD is chosen by the size of wavelet coefficients found from a wavelet 
decomposition of the numerical solution at a given time. One chooses a threshold 
with which to measure the coefficient size. That is, if the threshold is set to .001 
then the grid is refined at a given location and scale whenever the wavelet coefficients 
at that location and scale are larger in magnitude than .001. If the threshold is set 
to 0 then one gets finite difference on an evenly-spaced grid at the finest scale. The 
question then becomes, what is the relationship between this threshold value and the 
solution achieved by WOFD. As of now, a theoretical relationship does not exist but 
a numerical relationship does. That is, if the threshold is set to T = le~ p then one 
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can expect that both the h and /<*, difference between WOFD and finite difference 
at the finest scale will be a constant times this threshold, say kT where k < 10. For 
example, for a simulation with periodic boundary conditions, viscosity set at .01, and 
the final time set to 2, an h difference of 3.27e -4 and an difference of 2.31e -3 
were found for a threshold value of T = le~ 3 . This relationship was typical of all 
simulations which were run. 
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7 Stability of WOFD 


The discussion of the stability of WOFD will be given in terms of the eigenvalues of 
the differentiation matrix. 

7.1 Eigenvalues of Differentiation Matrix 

The differentiation matrix produced by this scheme will have a full set of eigenvectors. 
We can, therefore, look at instability through the magnitude of the eigenvalues. 

Recall that the WOFD method can produce essentially a completely arbitrary 
grid. The differentiation matrix can, therefore, take on an unlimited number of forms. 
For this reason, an analytical approach is not within reach. Therefore, experimental 
results which give the magnitude of the eigenvalues of the differentiation after each 
grid update will be given. On the following pages the eigenvalues will be given for 
the matrix, 


M = I + DA t + 1/2D 2 (A t) 2 + l/6D 3 (Ai) 3 + 1/24D 4 (A t)\ (48) 

which corresponds to WOFD being applied to the linear equation 


— 'U'X') 

with 4th-order Runge-Kutta time discretization. The grid is the grid that is chosen 
for the nonlinear Burgers’ equation. It is seen that the magnitude of the eigenvalues 
do sometimes exceed 1, but they rarely exceed 1 by very much. That is, for the 
periodic case, considering the maximum eigenvalue for the 4th-order RK for every 
grid encountered, the maximum eigenvalue magnitude for the entire run up to time 2 
is 1.0004. This eigenvalue is close enough to 1 in magnitude not to excite instability. 
In fact, the data would have to have a large component in the direction of the corre- 
sponding eigenvector and one would have to iterate 100 times to get amplification of 
4% in the direction of this eigenvector. 
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8 Efficiency of WOFD 

In a word, the efficiency of the WOFD method depends primarily on the data com- 
pression ratio. That is, we choose some finest grid which captures all the details of 
our solution throughout the entire run. Let’s say that this finest grid has N degrees- 
of-freedom. Now choose how close, in | • | 2 , you desire your WOFD solution to be 
to the solution on the finest grid. Choosing this ‘closeness parameter’ dictates the 
data compression ratio. Let’s say that the WOFD method needs only N 0 degrees- 
of-freedom to satisfy this criterion. Then, the amount of work done is, roughly, 
times the amount of work done to get the solution on the finest grid. 

8.1 Work Involved for Grid Update 

The grid update requires a number of steps. I will give a worst case scenario in 
estimating the number of operations. 

A grid update requires order N multiplies where N is the number of degrees-of- 
freedom in the finest scale. The constant that is multiplied times N is reasonably 
large, and the following will show where the operations are used: 

1. One must reconstruct the function on the finest grid. This requires about lOiV 
multiplies. 

2. Next, one must perform a wavelet decomposition. For a Daubechies 4 wavelet 
decomposition the filters are length 4. Therefore, the first decomposition re- 
quires 4iV multiplies. Likewise, the second decomposition requires 2 N multi- 
plies. The number of decompositions will determine the number of multiplies, 
but let’s say that this step, also, requires about lOiV multiplies. 

3. Choosing a grid from a wavelet decomposition does not require many operations, 
but it does need a number of ‘IF-THEN’ statements. There is roughly 1 ‘IF- 
THEN’ statement for each degree-of-freedom. Let me, once again, overestimate 
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the cost of each ‘IF-THEN’ and also compensate for a few other operations that 
are necessary by saying this requires roughly lOiV operations. 

4. Once a grid is chosen the new differentiation matrix is found. This requires 
about 40iVo multiplies in a worst case scenario. The second derivative filter can 
be found from the first derivative filter by a convolution of each filter. This 
requires about 25iVo multiplies. Here, N 0 is the number grid points used in the 
compressed scenario. N 0 is some fraction of N. 


The total number of multiplies is the sum of all the above multiples. That is, 
about 30 N + 65iVo. These numbers are rough, and we might as well round up to 
be safe and say, 50 N + lOOiVo multiplies are needed to define the grid and build 
a new differentiation matrix. This is reasonably expensive, but the update can be 
done rarely during a run. Compare this to finite difference on the finest scale. The 
filters for 1-st and 2-nd order differentiation are length 5. Each step of Runge-Kutta 
requires, therefore, at least lOiV multiplies. If we are using a fourth order RK then 
we have at least 40iV multiplies for each time step. It is, therefore, fair to say that 
the grid update step requires about the same amount of work as one time step taken 
using the full grid. 


8.2 Work Saved with Larger Time Step 


All the numerical scenarios use explicit time stepping. For Burgers’ equation with 
viscosity e the time step is set to, 


At 



1 


where Ax is the minimum spacing of the grid produce by WOFD. At the beginning of 
any simulation if the initial condition is smooth, as measured by a wavelet decomposi- 
tion, then the minimum Ax produced by WOFD is much larger than the Ax used on 
the finest grid. This allows a much larger time step without introducing large errors. 
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The larger time step means that far fewer total time steps will need to be taken to 
arrive at the final time. Fewer time steps gives a significant savings in terms of total 
operations performed. 
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9 WOFD in Higher Dimensions 


The effectiveness of WOFD has been illustrated in 1 dimension. The natural follow- 
up question would concern the effectiveness of WOFD in higher dimensions. At this 
point let us break WOFD into two parts: the first part is the grid definition, and the 
second part is the differencing on this new grid. 

Grid definition falls within the realm of local spectral analysis. That is, one is 
interested in the spectrum locally. Local high frequency data requires a grid density 
sufficient to resolve the highest frequency, whereas local low frequency data can be re- 
solved with a relatively coarse grid. The wavelet structure provides a very convenient 
mechanism to perform this nested group of short-time Fourier transforms. For higher 
dimensions, say 2 dimensions, one need only choose a coordinate system for the space 
and simply perform the wavelet filtering throughout all of the data and parallel to 
each axis. The most common question at this point concerns the effectiveness of this 
method of grid selection when the data is composed of structure which is at a 45 
degree angle to the axes. The simple answer is that it works well since all structures 
within the data can be projected onto the orthogonal coordinate system which spans 
the space. If, however, one is not satisfied with the grid given in this situation then 
the parameter which adjusts the sensitivity of grid selection can be adjusted. With 
this sensitivity adjustment, one will always fold a suitable grid. Included here are 
three sets of data which will illustrate the effective grid selection. The first function 
is a discontinuity at a 45 degree angle to the axes, see Figure (8), and the grid chosen 
for this data, see (9). The second set of data corresponds, roughly, to the inner and 
outer flow near a boundary, 

f(x,y) = l-e~2~, (49) 

see Figure (10), and the grid chosen in this case, see Figure (11). The third set of data 
is numerically-generated pressure from a turbulent flow. A contour diagram is given 
in Figure (12) and the grid chosen is given in Figure (13). In all cases the wavelet is 
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the D 4 wavelet. 

Differencing on a grid chosen by WOFD will depend on the application. The D 4 is 
the optimal wavelet if one is using a central optimal 4th-order finite difference method 
in the sense of Section 3. But one is not confined to matching the wavelet precisely 
to the differencing method. The only recommendation is that the stencil of the finite 
difference method be of roughly the same size as the length of the wavelet filters. 
This will insure that in the grid selection one is not filtering data which is too far 
outside of the support of the differencing stencil. Also, depending on the application 
one might need a reliable mechanism for choosing locally rectangular grids. In this 
case one would choose the finest grid produced by WOFD in each rectangular region. 

One other issue is the use of WOFD to refine beyond the ‘finest grid’. For the 1 
dimension case considered in this paper there was always an underlying finest grid. 
This was a theoretical convenience and not a necessity or even a desirable feature. 
Based on the energy or magnitude of the smallest scale wavelet coefficients one can 
refine to any desirable grid density, see [10]. 
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Function which is decomposed by Wavelets 



Figure 8: A discontinuity at a 45 degree angle to the axes. 
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Function which is decomposed by Wavelets 



Figurfr 10: A snapshot of the inner and outer flow near a boundary. 
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Grid Chosen by Wavelets 



Figure 11: The grid chosen by the wavelet decomposition of the flow near a boundary. 
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0 





Figure 13: The grid chosen by the wavelet decomposition of the pressure from tur- 


bulent data. 
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10 Conclusion 


The WOFD method is essentially the same as a Daubechies-based wavelet method. 
However, two serious problems encountered with wavelet numerical methods are over- 
come. That is, boundary conditions can now be imposed in exactly the same manner 
they are imposed for finite difference methods. Furthermore, there is no longer a 
problem with nonlinear terms since we are now working with the point values of the 
function in the physical space. WOFD will approximate a conservative numerical 
method as long as one is working with a conservative numerical method on the finest 
uniform grid. This approximation can be as fine as the user wishes. 

In this paper the WOFD method has been explored for the case of a 5-point 4- 
th order finite difference operator on an arbitrary grid. We can, however, use this 
method with higher order schemes by using the filters associated with the higher 
order Daubechies wavelets to define the grid. The only suggestion is that the stencil 
size of the numerical scheme be of roughly the same size as the length of the filters. 

Higher dimensional applications of WOFD remain to be explored. As mentioned, 
WOFD can be broken into two parts, the grid selection and the differencing. The 
grid selection step requires the wavelet analysis to detect the local oscillation content 
of the data. Based on this grid, one can choose a number of ways to apply finite 
differencing. It has been shown here that grid selection is quite effective for ‘difficult’ 
data sets in two dimensions. Future research will combine this 2-dimensional grid 
selection with appropriate differencing. 
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